Transformación

Sabemos que una de las formas de estabilizar la variancia en nuestra serie es realizando alguna transformación como lo son la logarítmica o cuadrática. Se probaron algunas transformaciones en la serie pero ninguna parecia tener una implicancia positiva en la estabilizacion de la variancia, por ello se decidió presentar la transformacion recomendad por el gráfico de Box y Cox solo a modo de ejemplo

boxcox(lm(data = emae_mensual_ts, emae_original ~ 1), 
       plotit = T, 
       lambda = seq(-3, 8, by = 0.5))

lambda <- 3
# Aplicar la transformación
emae_mensual_ts <- emae_mensual_ts %>%
  mutate(emae_transformed = (emae_original^lambda - 1) / lambda)

ggplot(emae_mensual_ts, aes(x = year_factor, y = emae_transformed)) +
  geom_boxplot(fill="cyan4") +
  labs(title = "Boxplot de Valores anuale", x = "Year", y = 'Indice 2004 = 100') +
  theme_minimal() # no tiene variancia constante

Como mencionamos anteriormente, esta transformación solo parece cambiar la escala pero no la forma de la variancia anual, por lo que continuaremos nuestro trabajo con la serie original.

Diferenciaciones

Al analizar la serie temporal del EMAE, observamos un patrón de estacionalidad recurrente, donde las fluctuaciones anuales se mantienen consistentemente a lo largo del tiempo. Para abordar este comportamiento y trabajar con una serie estacionaria (requisito fundamental para los modelos ARIMA), decidimos aplicar una diferenciación estacional. Este proceso nos permite eliminar los componentes estacionales y centrar nuestro análisis en las variaciones más relevantes de la tendencia y los ciclos.

emae_mensual_ts %>% 
  gg_tsdisplay(difference(emae_original, 12), 
               plot_type='partial',
               lag=36)

Al analizar los gráficos de la Función de Autocorrelación (FAC) y la Función de Autocorrelación Parcial (FACP), identificamos características que sugieren un modelo autorregresivo (AR) en la parte regular de la serie. En particular:

• La FAC muestra una caída exponencial y un patrón sinusoidal, característico de un modelo AR.

• La FACP resalta significancia en los rezagos 2 y 3, lo que podría apuntar a un modelo \(AR(3)\). Sin embargo, optamos por el modelo más sencillo, \(AR(1)\), para mantener la parsimonia.

Para la parte estacional, observamos picos significativos en los rezagos 12 y 24, lo que indica una estacionalidad anual bien definida. Luego de aplicar una diferenciación estacional, consideramos los siguientes modelos como candidatos:

• \(ARIMA(1,1,0)\): Un modelo simple que captura la dinámica principal.

• \(ARIMA(2,1,0)\): Una extensión que podría reflejar interacciones adicionales entre rezagos.

Finalmente, la selección del modelo se basará en criterios de ajuste como AIC, BIC y análisis de residuos para garantizar un balance la bondad de ajuste y lo parsimoñoso del modelo.

Ajuste de la Serie

A partir del análisis previo de los gráficos de la FAC y la FACP, identificamos posibles configuraciones para los modelos ARIMA que se ajustan tanto a la parte regular como a la estacional de la serie. Con base en estos hallazgos, y para garantizar un análisis robusto, procederemos a ajustar los modelos identificados de dos maneras:

• Modelos identificados manualmente: Utilizaremos las configuraciones sugeridas (por ejemplo, \(ARIMA(1,1,0)\) y \(ARIMA(2,1,0)\)) como punto de partida para verificar su desempeño en el ajuste de la serie.

• Modelos automáticos: Exploraremos opciones generadas automáticamente mediante métodos de selección óptimos basados en criterios como AIC y BIC, que evalúan el balance entre la complejidad del modelo y su capacidad de predicción.

De esta manera, compararemos los resultados de ambos enfoques para determinar cuál ofrece el mejor ajuste y predicción, siempre considerando la parsimonia y la calidad de los residuos.

# Definir nombres de archivo para guardar los modelos
fit_original_file <- "fitted/fit_original.rds"
fit_transformed_file <- "fitted/fit_transformed.rds"

# Ajuste y guardado de fit_original
NOTFIT = TRUE
if (NOTFIT) {
  fit_original <- readRDS(fit_original_file)
  message("Modelo fit_original cargado desde \fitted.")
} else {
  fit_original <- emae_mensual_ts %>% 
    model(
      arima11 = ARIMA(emae_original ~ pdq(1,0,0) + PDQ(1,1,0)),
      arima11_2 = ARIMA(emae_original ~ pdq(1,0,0) + PDQ(2,1,0)),
      auto = ARIMA(emae_original, stepwise = FALSE, approx = FALSE),
      auto_2 = ARIMA(emae_original, stepwise = TRUE, approx = TRUE)
    )
  saveRDS(fit_original, fit_original_file)
  message("Modelo fit_original ajustado y guardado en \fitted.")
}


df_fit <- as.data.frame(fit_original %>% 
                           pivot_longer(everything(),
                                        names_to = "Model name",
                                        values_to ="Orders"))  
kable(df_fit,
      format = "html", 
      table.attr = "class='table table-striped'") %>%
  kable_styling() %>%
  row_spec(0, bold = T, background = "cyan4", color = "white")

glance_summary <- glance(fit_original) %>%
  select(.model, sigma2, log_lik, AIC, AICc, BIC)

kable(as.data.frame(glance_summary), 
      format = "html", 
      table.attr = "class='table table-striped'") %>%
  kable_styling() %>%
  row_spec(0, bold = T, background = "cyan4", color = "white")

De la tabla presentada, obtenemos las siguientes observaciones clave:

1. El mejor modelo identificado es el ajustado automáticamente, \(SARIMA(0,1,2)(2,1,2)_{12}\). Este modelo utiliza una diferenciación y un comportamiento \(MA(2)\) en la parte regular y una diferenciación y un comportamiento \(ARMA(2,2)\) en la parte estacional.

2. El modelo \(SARIMA(1,0,0)(2,1,0)_{(12)}\) presenta un desempeño competitivo y será considerado como una alternativa más sencilla para comparación.

A partir de estas conclusiones, trabajaremos con ambos modelos (\(SARIMA(0,1,2)(2,1,2)_{12}\). y \(SARIMA(1,0,0)(2,1,0)_{(12)}\)) en las siguientes etapas del análisis, evaluando su desempeño en términos de ajuste y predicción.

params_df <- as.data.frame(fit_original %>%
                              tidy() %>%  
                              select(.model, term, estimate)%>%
                              filter(.model %in% c("auto", "arima11_2"))
                          )

kable(as.data.frame(params_df), 
      format = "html", 
      table.attr = "class='table table-striped'") %>%
  kable_styling() %>%
  row_spec(0, bold = T, background = "cyan4", color = "white")

fitted_values <- fit_original %>%
  augment() %>%  # Obtener valores ajustados
  as_tibble() %>% # Convertir a tibble para ggplot
  mutate(Modelo = .model)

fitted_values <- fitted_values %>%
  mutate(periodo_date = parse_date_time(periodo, orders = "Y b")) # Convierte fecha

fitted_values_filtered <- fitted_values %>%
  filter(
    periodo_date >= ymd("2019-01-01") & periodo_date <= ymd("2024-12-31"),
    .model %in% c("auto", "arima11_2") # Filtrar modelos
  )

# Graficar los valores filtrados
g <- ggplot(data = fitted_values_filtered, aes(x = periodo_date)) +
  geom_line(aes(y = emae_original, color = "Serie Original"), size = 1) +
  geom_line(aes(y = .fitted, color = .model), linetype = "dashed", size = 0.8) +
  labs(
    title = "Serie Original y Estimaciones de los Modelos (2019-2024)",
    x = "Tiempo",
    y = "Índice 2004 = 100",
    color = "Leyenda"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("Serie Original" = "cyan4", "arima11_2" = "#009933", "auto" = "#9933CC"))
ggplotly(g)

Esperanza y Variancia de los residuos

Comenzaremos con un gráfico de Residuos vs Timepo para comprobar esperanza y variancia constante

# Extraer los residuos del modelo "auto"
residuos1 <- residuals(fit_original, type = "response") %>%
  filter(.model == "auto")  #  modelo "auto"
residuos2 <- residuals(fit_original, type = "response") %>%
  filter(.model == "arima11_2")  #arima11_2

# Graficar los residuos vs tiempo
g1 <- ggplot(residuos1, aes(x = periodo, y = .resid)) +
  geom_point(color = "#9933CC") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Residuos Modelo Automatico",
    x = "",
    y = "Residuos"
  ) +
  theme_minimal()
g2 <- ggplot(residuos2, aes(x = periodo, y = .resid)) +
  geom_point(color = "#009933") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Residuos Modelo Manual",
    x = "",
    y = "Residuos"
  ) +
  theme_minimal()
g1+g2

Se observan valores atípicos evidentes en marzo y abril de 2020, correspondientes al inicio de la pandemia de COVID-19. Excluyendo estos puntos, los residuos del modelo automático parecen mostrar una varianza constante. En el caso del modelo identificado “manualmente,” se aprecia un ligero incremento en la varianza que, sin embargo, podríamos despreciar. Esta estabilidad en los residuos nos habla de dosbuenos ajustes.

Independencia de los Residuos

Continuamos nuestro análisis con uno el más importante de los supuestos, la incorrelación de los residuos. Para evaluarlo, observamos los gráficos de autocorrelación (FAC) de los residuos de ambos modelos.

# Gráfico de autocorrelación para los residuos
# acf(residuos1$.resid, main = "Autocorrelación de los residuos", lag.max = 48)
# acf(residuos2$.resid, main = "Autocorrelación de los residuos", lag.max = 48)

acf1_data <- acf(residuos1$.resid, plot = FALSE)
acf2_data <- acf(residuos2$.resid, plot = FALSE)
acf_df <- data.frame(
  Lag1 = acf1_data$lag,
  ACF1 = acf1_data$acf,
  Lag2 = acf2_data$lag,
  ACF2 = acf2_data$acf
)

# Graficar autocorrelación con ggplot
r1 <- ggplot(acf_df, aes(x = Lag1, y = ACF1)) +
  geom_bar(stat = "identity", fill = "#9933CC", alpha = 0.7, width = 0.2) +
  geom_hline(yintercept = 0, color = "black") +
  geom_hline(yintercept = c(-1.96 / sqrt(length(residuos1$.resid)),
                            1.96 / sqrt(length(residuos1$.resid))),
             linetype = "dashed", color = "red") +
  labs(title = "Autocorrelación Modelo Automatico", 
       x = "Lags",
       y = "Autocorrelación de residuos") +
  theme_minimal()
r2 <- ggplot(acf_df, aes(x = Lag2, y = ACF2)) +
  geom_bar(stat = "identity", fill = "#009933", alpha = 0.7, width = 0.2) +
  geom_hline(yintercept = 0, color = "black") +
  geom_hline(yintercept = c(-1.96 / sqrt(length(residuos2$.resid)),
                            1.96 / sqrt(length(residuos2$.resid))),
             linetype = "dashed", color = "red") +
  labs(title = "Autocorrelación Modelo Manual", 
       x = "Lags",
       y = "Autocorrelación de residuos") +
  theme_minimal()
r1+r2

#TODO: indica errores correlacionas para la parte estacional?

En el caso del modelo automático, la FAC muestra un comportamiento consistente con un ruido blanco, con la excepción de un bastón significativo en un rezago mayor al 20. Esto sugiere que, en general, los residuos están incorrelados.

Por otro lado, la FAC de los residuos del modelo identificado manualmente revela varios rezagos significativos antes del lag 15, lo que indica que estos residuos no se comportan completamente como un ruido blanco. Este hallazgo podría reflejar deficiencias en el ajuste de este modelo, especialmente al momento de realizar predicciones.

#TODO: test test de ljung box para los rezagos

Normalidad de los Residuos

Finalmente, observamos un histograma y un Q-Q plot de los residuos para evaluar su proximidad a una distribución normal.

# Histograma de los residuos
h1 <- ggplot(residuos1, aes(x = .resid)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "#9933CC", alpha = 0.7) +
  geom_density(color = "red", size = 1) +
  labs(
    title = "Histograma de modelo automatico",
    x = "Residuos",
    y = "Densidad"
  ) +
  theme_minimal()

# QQ-Plot normal
qq1 <- ggplot(residuos1, aes(sample = .resid)) +
  stat_qq(color = "#9933CC") +
  stat_qq_line(color = "red", linetype = "dashed") +
  labs(
    title = "QQ-Plot de los residuos",
    x = "Cuantiles teoricos",
    y = "Cuantiles de los residuos"
  ) +
  theme_minimal()

h2 <- ggplot(residuos2, aes(x = .resid)) +
  geom_histogram(aes(y = ..density..), bins = 30, fill = "#009933", alpha = 0.7) +
  geom_density(color = "red", size = 1) +
  labs(
    title = "Histograma de modelo manual",
    x = "Residuos",
    y = "Densidad"
  ) +
  theme_minimal()

# QQ-Plot normal
qq2 <- ggplot(residuos2, aes(sample = .resid)) +
  stat_qq(color = "#009933") +
  stat_qq_line(color = "red", linetype = "dashed") +
  labs(
    title = "QQ-Plot de los residuos",
    x = "Cuantiles teoricos",
    y = "Cuantiles de los residuos"
  ) +
  theme_minimal()

(h1+h2) / (qq1+qq2)

Si bien identificamos la presencia de outliers (Marzo y Abril de 2020), los cuales desvían la distribución de los residuos de la normalidad, esto no es motivo de preocupación ya que el supuesto de normalidad no es estrictamente necesario para la calidad del modelo. Utilizamos este supuesto para construir intervalos de confianza en nuestras prediccione.

Medidas de error

Para evaluar y comparar las predicciones de nuestros modelos utilizaremos el MAPE (Mean Absolute Percentage Error). Este representa el promedio del error absoluto en términos porcentuales respecto a los valores reales, indica, en promedio, cuánto se desvió porcentualmente el modelo de los valores reales.

La fórmula del MAPE es:

\[ MAPE = \frac{1}{M} \sum_{l=1}^M \left| \frac{e_l}{Z_{n+l}} \right| \times 100 \] Donde:

• \(e_l =Z_{n+l} - \hat{z_n(l)}\) error de la predicción.

Un MAPE bajo indica un buen rendimiento, pero esta métrica puede ser sensible a valores reales cercanos a cero, lo que puede inflar el porcentaje de error, por esto la acompañaremos con el MSE (Mean Square Error).

Utilizamos las funciones MAPE y MSE de la libreria MLmetrics

mape_auto <- round(MAPE(y_pred = values_to_predict, y_true = f1$.mean)*100,2)
mse_auto <- MSE(y_pred = values_to_predict, y_true = f1$.mean)
mape_manual <- round(MAPE(y_pred = values_to_predict, y_true = f2$.mean)*100,2)
mse_manual <- MSE(y_pred = values_to_predict, y_true = f2$.mean)

metricas <- data.frame(
  Modelo = c("Automático", "Manual"),
  MAPE = c(mape_auto, mape_manual),
  MSE = c(mse_auto, mse_manual)
)

# Mostrar como tabla con kable
kable(metricas, 
      col.names = c("Modelo", "MAPE", "MSE"), 
      caption = "Metricas para comparar modelos", 
      format = "html", 
      table.attr = "class='table table-striped'") %>%
  kable_styling() %>%
  row_spec(0, bold = T) %>%
  row_spec(1, bold = T, background = "#9933CC", color = "white") %>%
  row_spec(2, bold = T, background = "#009933", color = "white")

Obervando estas metricas podemos decir que ambos modelos muestran un excelente ajuste a los datos originales. En promedio, los modelos se equivocan en apenas por encima del 2% al realizar predicciones, lo cual nos habla de una gran precisión predictiva. Podemos decir que cualquiermodelo es confiables para la predicción, siendo el modelo identificado “a ojo” más preciso y más parsimonioso.

Es importante destacar que, si bien estamos prediciendo valores desconocidos para la serie, estos no son desconocidos para nosotros. Podríamos haber entrenado el modelo con datos hasta Agosto de 2024 y realizar predicciones para los valores futuros del EMAE.